load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

location="/Users/ivana/WORK_BSC/final_figures_allmodels/ecmwf/"


fc    = addfile(location+"g250hpa_era5_19792004djf.nc","r")
fw    = addfile(location+"g250hpa_era5_20092019djf.nc", "r")

;*******************************************;
;Read in variables
;*******************************************

lat=fc->latitude
lon=fc->longitude

zc=short2flt(fc->z(:,0,:,:))
zw=short2flt(fw->z(:,0,:,:))

;***********************************************
; average over the simulations/years
;************************************************

tcs=dim_avg(zc(latitude|:,longitude|:,time|:))
tws=dim_avg(zw(latitude|:,longitude|:,time|:))

;************************
; calculate the anomaly:
;************************

  finfo_ano=tws
  finfo_ano=(tws - tcs)/9.81

  finfo_ano!0 = "lat"
  finfo_ano&lat = lat

  finfo_ano!1 = "lon"
  finfo_ano&lon =lon

;***********************************
; perform t-test
;***********************************

  xtmp = zc(latitude|:,longitude|:,time|:)     ; reorder but do it only once [temporary]
  ytmp = zw(latitude|:,longitude|:,time|:)

  xAve = dim_avg (xtmp)              ; calculate means at each grid point
  yAve = dim_avg (ytmp)
  xVar = dim_variance (xtmp)         ; calculate variances
  yVar = dim_variance (ytmp)
  sigr = 0.1                        ; critical sig lvl for r
  xEqv = equiv_sample_size (xtmp, sigr,0)
  yEqv = equiv_sample_size (ytmp, sigr,0)
;
   xN=xEqv
   yN=yEqv
   iflag= True                         ; population variance not similar
 ;iflag= False
  prob = ttest(xAve,xVar,xN, yAve,yVar,yN, iflag, False)
  prob!0 = "lat"
  prob!1 = "lon"
  prob&lat =lat
  prob&lon =lon

;-------------------------
; making the plots
;-------------------------
wks = gsn_open_wks("pdf","zg250hpa_ano_djf_era5_20092019min19792004_ttest90")

gsn_define_colormap(wks,"BlueDarkOrange18")

plot = new(1,graphic)                          ; create a plot array

  res          = True
  res@gsnDraw  = False                          ; don't draw
  res@gsnFrame = False                          ; don't advance frame, these 2 are going to be set with panel plot - gsn_panel!
  res@cnInfoLabelOn = False                     ; turn off cn info label
  res@cnLinesOn           = False
  res@cnFillOn            = True          ; turn on color
  res@lbLabelBarOn        = True           ; turn off individual cb's
 
  res@mpPerimOn = False
  res@cnLevelSelectionMode =  "ExplicitLevels"

 res@cnLineLabelsOn=False
 res@mpCenterLonF     = 300. 
 res@mpGeophysicalLineThicknessF = 2.

 res@gsnPolar   = "NH"
 res@mpMinLatF            = 30

  res1          = True
  res1@gsnDraw  = False                          ; don't draw
  res1@gsnFrame = False                          ; don't advance frame, these 2 are going to be set with panel plot - gsn_panel!
  res1@cnInfoLabelOn = False                     ; turn off cn info label
  res1@cnLinesOn           = False; True
  res1@cnFillOn            = True          ; turn on color
  res1@lbLabelBarOn        = False           ; turn off individual cb's

  res1@cnMonoFillPattern = False ;True
  res1@cnFillPatterns      = (/-1,8,-1/) ;17, 8, 13
  res1@cnLevelSelectionMode =  "ExplicitLevels"
  res1@cnLevels =(/0.0, 0.1/)  ; ! 95%, 90% level

  res1@cnMonoFillColor     = True         ; Use same fill color
  res1@cnFillColor         = (/"black"/)

  res1@cnLineLabelsOn=False
  res1@cnMonoLineColor= True
  res1@cnLineColor      = (/"grey36"/) ; (/"white"/)

;res@mpGeophysicalLineThicknessF = 1.5

res1@gsnCenterString = ""
res1@gsnLeftString = ""
res1@gsnRightString = ""

res@gsnLeftString = "DJF"
res@gsnLeftString = ""
res@gsnCenterString = ""

res@cnFillColors = (/4,5,10,0,0,12,13,14,15,16,17,18/)
res@cnLevels =(/-20.,-10.,-5.,0,5.,10.,20.,30.,40.,50.,60./)


plot(0) = gsn_csm_contour_map(wks,finfo_ano(:,:),res)
plot2 = gsn_csm_contour(wks,gsn_add_cyclic_point(prob),res1)
overlay(plot(0),plot2)

;________________________
; creating the panel:
;________________________
 resP            = True

resP@gsnPanelLabelBar    = False;True                ; add common colorbar
resP@lbLabelFontHeightF  = 0.007               ; make labels smaller

resP@gsnPanelBottom   = 0.05                   ; add space at bottom
;resP@gsnPanelFigureStrings= (/"a)","b)","c)"/) ; add strings to panel
res@txFontHeightF     = .24


resP@gsnMaximize    = True                ; maximize plot
gsn_panel(wks,plot,(/1,1/),resP)


